<?php

/**
 * 	@package JAMA
 *
 * 	For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
 * 	orthogonal matrix Q and an n-by-n upper triangular matrix R so that
 * 	A = Q*R.
 *
 * 	The QR decompostion always exists, even if the matrix does not have
 * 	full rank, so the constructor will never fail.  The primary use of the
 * 	QR decomposition is in the least squares solution of nonsquare systems
 * 	of simultaneous linear equations.  This will fail if isFullRank()
 * 	returns false.
 *
 * 	@author  Paul Meagher
 * 	@license PHP v3.0
 * 	@version 1.1
 */
class PHPExcel_Shared_JAMA_QRDecomposition
{

    const MatrixRankException = "Can only perform operation on full-rank matrix.";

    /**
     * 	Array for internal storage of decomposition.
     * 	@var array
     */
    private $QR = array();

    /**
     * 	Row dimension.
     * 	@var integer
     */
    private $m;

    /**
     * 	Column dimension.
     * 	@var integer
     */
    private $n;

    /**
     * 	Array for internal storage of diagonal of R.
     * 	@var  array
     */
    private $Rdiag = array();

    /**
     * 	QR Decomposition computed by Householder reflections.
     *
     * 	@param matrix $A Rectangular matrix
     * 	@return Structure to access R and the Householder vectors and compute Q.
     */
    public function __construct($A)
    {
        if ($A instanceof PHPExcel_Shared_JAMA_Matrix)
        {
            // Initialize.
            $this->QR = $A->getArrayCopy();
            $this->m  = $A->getRowDimension();
            $this->n  = $A->getColumnDimension();
            // Main loop.
            for ($k = 0; $k < $this->n; ++$k)
            {
                // Compute 2-norm of k-th column without under/overflow.
                $nrm = 0.0;
                for ($i = $k; $i < $this->m; ++$i)
                {
                    $nrm = hypo($nrm, $this->QR[$i][$k]);
                }
                if ($nrm != 0.0)
                {
                    // Form k-th Householder vector.
                    if ($this->QR[$k][$k] < 0)
                    {
                        $nrm = -$nrm;
                    }
                    for ($i = $k; $i < $this->m; ++$i)
                    {
                        $this->QR[$i][$k] /= $nrm;
                    }
                    $this->QR[$k][$k] += 1.0;
                    // Apply transformation to remaining columns.
                    for ($j = $k + 1; $j < $this->n; ++$j)
                    {
                        $s = 0.0;
                        for ($i = $k; $i < $this->m; ++$i)
                        {
                            $s += $this->QR[$i][$k] * $this->QR[$i][$j];
                        }
                        $s = -$s / $this->QR[$k][$k];
                        for ($i = $k; $i < $this->m; ++$i)
                        {
                            $this->QR[$i][$j] += $s * $this->QR[$i][$k];
                        }
                    }
                }
                $this->Rdiag[$k] = -$nrm;
            }
        }
        else
        {
            throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);
        }
    }

//	function __construct()

    /**
     * 	Is the matrix full rank?
     *
     * 	@return boolean true if R, and hence A, has full rank, else false.
     */
    public function isFullRank()
    {
        for ($j = 0; $j < $this->n; ++$j)
        {
            if ($this->Rdiag[$j] == 0)
            {
                return false;
            }
        }
        return true;
    }

//	function isFullRank()

    /**
     * 	Return the Householder vectors
     *
     * 	@return Matrix Lower trapezoidal matrix whose columns define the reflections
     */
    public function getH()
    {
        for ($i = 0; $i < $this->m; ++$i)
        {
            for ($j = 0; $j < $this->n; ++$j)
            {
                if ($i >= $j)
                {
                    $H[$i][$j] = $this->QR[$i][$j];
                }
                else
                {
                    $H[$i][$j] = 0.0;
                }
            }
        }
        return new PHPExcel_Shared_JAMA_Matrix($H);
    }

//	function getH()

    /**
     * 	Return the upper triangular factor
     *
     * 	@return Matrix upper triangular factor
     */
    public function getR()
    {
        for ($i = 0; $i < $this->n; ++$i)
        {
            for ($j = 0; $j < $this->n; ++$j)
            {
                if ($i < $j)
                {
                    $R[$i][$j] = $this->QR[$i][$j];
                }
                elseif ($i == $j)
                {
                    $R[$i][$j] = $this->Rdiag[$i];
                }
                else
                {
                    $R[$i][$j] = 0.0;
                }
            }
        }
        return new PHPExcel_Shared_JAMA_Matrix($R);
    }

//	function getR()

    /**
     * 	Generate and return the (economy-sized) orthogonal factor
     *
     * 	@return Matrix orthogonal factor
     */
    public function getQ()
    {
        for ($k = $this->n - 1; $k >= 0; --$k)
        {
            for ($i = 0; $i < $this->m; ++$i)
            {
                $Q[$i][$k] = 0.0;
            }
            $Q[$k][$k] = 1.0;
            for ($j = $k; $j < $this->n; ++$j)
            {
                if ($this->QR[$k][$k] != 0)
                {
                    $s = 0.0;
                    for ($i = $k; $i < $this->m; ++$i)
                    {
                        $s += $this->QR[$i][$k] * $Q[$i][$j];
                    }
                    $s = -$s / $this->QR[$k][$k];
                    for ($i = $k; $i < $this->m; ++$i)
                    {
                        $Q[$i][$j] += $s * $this->QR[$i][$k];
                    }
                }
            }
        }
        /*
          for($i = 0; $i < count($Q); ++$i) {
          for($j = 0; $j < count($Q); ++$j) {
          if(! isset($Q[$i][$j]) ) {
          $Q[$i][$j] = 0;
          }
          }
          }
         */
        return new PHPExcel_Shared_JAMA_Matrix($Q);
    }

//	function getQ()

    /**
     * 	Least squares solution of A*X = B
     *
     * 	@param Matrix $B A Matrix with as many rows as A and any number of columns.
     * 	@return Matrix Matrix that minimizes the two norm of Q*R*X-B.
     */
    public function solve($B)
    {
        if ($B->getRowDimension() == $this->m)
        {
            if ($this->isFullRank())
            {
                // Copy right hand side
                $nx = $B->getColumnDimension();
                $X  = $B->getArrayCopy();
                // Compute Y = transpose(Q)*B
                for ($k = 0; $k < $this->n; ++$k)
                {
                    for ($j = 0; $j < $nx; ++$j)
                    {
                        $s = 0.0;
                        for ($i = $k; $i < $this->m; ++$i)
                        {
                            $s += $this->QR[$i][$k] * $X[$i][$j];
                        }
                        $s = -$s / $this->QR[$k][$k];
                        for ($i = $k; $i < $this->m; ++$i)
                        {
                            $X[$i][$j] += $s * $this->QR[$i][$k];
                        }
                    }
                }
                // Solve R*X = Y;
                for ($k = $this->n - 1; $k >= 0; --$k)
                {
                    for ($j = 0; $j < $nx; ++$j)
                    {
                        $X[$k][$j] /= $this->Rdiag[$k];
                    }
                    for ($i = 0; $i < $k; ++$i)
                    {
                        for ($j = 0; $j < $nx; ++$j)
                        {
                            $X[$i][$j] -= $X[$k][$j] * $this->QR[$i][$k];
                        }
                    }
                }
                $X = new PHPExcel_Shared_JAMA_Matrix($X);
                return ($X->getMatrix(0, $this->n - 1, 0, $nx));
            }
            else
            {
                throw new PHPExcel_Calculation_Exception(self::MatrixRankException);
            }
        }
        else
        {
            throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
        }
    }

//	function solve()
}

//	PHPExcel_Shared_JAMA_class QRDecomposition
